Quasiharmonic elastic constants corrected for deviatoric thermal 

stresses 



Pierre Carrier and Renata M. Wentzcovitch 
Minnesota Supercomputing Institute, 
Department of Chemical Engineering and Materials Science, 
University of Minnesota, Minneapolis, MN 55455 

Joao F. Justo 
Escola Politecnica, Universidade de Sao Paulo, 
CP 61548, CEP 05424-970, Sao Paulo, Brazil, 
and Chemical Engineering and Materials Science Department, 
University of Minnesota, MN 55455 
(Dated: August 25, 2008) 

Abstract 

The quasiharmonic approximation (QHA), in its simplest form also called the statically con- 
strained (SC) QHA, has been shown to be a straightforward method to compute thermoelastic 
properties of crystals. Recently we showed that for non-cubic solids SC-QHA calculations develop 
deviatoric thermal stresses at high temperatures. Relaxation of these stresses leads to a series 
of corrections to the free energy that may be taken to any desired order, up to self-consistency. 
Here we show how to correct the elastic constants obtained using the SC-QHA. We exemplify the 
procedure by correcting to first order the elastic constants of MgSi03-perovskite and MgSi03-post- 
perovskite, the major phases of the Earth's lower mantle. We show that this first order correction 
is quite satisfactory for obtaining the aggregated elastic averages of these minerals and their ve- 
locities in the lower mantle. This type of correction is also shown to be applicable to experimental 
measurements of elastic constants in situations where deviatoric stresses can develop, such as in 
diamond anvil cells. 

PACS numbers: 62.20.Dc, 65.40.-b, 91.35.-x, 91.60.Gf 
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I. INTRODUCTION 



The quasiharmonic approximation (QHA)W is a computationally efficient method for 
evaluating thermal properties of materials within the density functional theory (DFT) from 
low to temperatures above the Debye temperature. It provides high quality high pressure- 
high temperature materials properties^ 1 ^ 1 ^ in a continuous pressure-temperature (PT) 
domain in which anharmonic effects are negligible.- However, it has a not well recognized 
shortcoming: the non-hydrostatic nature of thermal stresses in non-isotropic structures. 
Broadly speaking, these calculations start by obtaining the static internal energy of fully 
relaxed DFT structures at various pressures. After computations of the vibrational density 
of states, the thermal energy contribution to the Helmholtz free energy is added. This latter 
contribution has anisotropic strain gradients and produces deviatoric stresses. This straight- 
forward procedure should be referred to as the statically constrained (SC) QHA. It has been 
used to compute the elastic constant tensor of isotropic^ and non-isotropic minerals-^ at 
high PT as well, even though pressure conditions were not precisely hydrostatic in the latter 
calculations. In general, relaxation of deviatoric stresses, irrespective of their origin, is essen- 
tial in both experiments^ 1 ^ and theory- for generating realistic and reproducible structural 
and elastic properties. 

Here we show how to correct the elastic constant tensor obtained using the SC-QHA. 
We exemplify the procedure by correcting to first order the elastic constants of MgSiCv 
perovskite (PV) and MgSi03-post-perovskite (PPV), the major phases of the Earth's lower 
mantle, for which elasticity data are essential to interpret seismic information of this region.— 
We show that this first order correction is quite satisfactory for obtaining the aggregated 
elastic averages of these minerals and their acoustic velocities in the PT range of the lower 
mantle. 

This article is organized as follow: we first discuss the equations used for numerically 
determining the elastic constant tensor within the SC-QHA. We then describe the proce- 
dure for correcting it to first order for deviatoric thermal stresses. We then evaluate these 
corrections to the previously reported elastic constant tensors of PV— and PPV. 6 
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II. ELASTICITY WITHIN AND BEYOND THE STATICALLY CONSTRAINED 
(SC) QHA 



The present procedure builds on a related procedure to correct structural parameters and 
equations of state of non-isotropic solids at high PTs. 9 The method introduced in Ref. can 
correct the SC crystal structure at V(P, T) to infinite order as long as the SC elastic constant 
tensor is simultaneously corrected. However, this is a very demanding computational pro- 
cedure and, fortunately, unnecessary. A first order correction to the crystal structure using 
SC elastic constant, appears to be sufficient. This conclusion was reached after examining 
the crystal structure of one of the most studied materials at high PT: MgSiOa-perovskite. 12 
This type of experimental data is quite limited and results on other materials with similarly 
complex crystal structures would be helpful strengthen this conclusion. 

According to the (SC) QHA the Helmholtz free energy is given by: 



F(V,T) 



E{v)+ j2 huj ^ v) 
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where kg and h are respectively Boltzmann's and Planck's constants. The first term, E(V), 
is the volume dependent static energy obtained after full structural relaxation under isotropic 
pressure, and ou(V) is the corresponding phonon spectrum. Both phonon spectrum and static 
energy are here determined using the DFT within the local density approximation (LDA),— 
but the methodology is general and applicable to any first principles method. Structural 
relaxations are performed using a variable cell shape (VCS) algorithm^ and phonon spectra 
are computed using the PWscf code^ as described in Ref. |l6|, based the linear response 
theory. The second term in Eq. ([T]) is the zero point energy, Fzp, such that the sum of 
the terms in the bracket is the energy at T=0 K. The last term in Eq. ([1]) is the thermal 
excitation energy, F t h (see Ref. O for details). 



Pressure, P, is obtained from F using the standard thermodynamics relation 

(2) 
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This procedure implicitly assumes that P remains isotropic at all temperatures, but this is 
only true for static calculations, where structures were optimized at target pressures. The 
two frequency dependent terms in Eq. ([1]), the zero-point energy and the thermal energy, 
contribute to P but their strain gradients are intrinsically anisotropic. This effect was 



recently quantified 9 by the computation of deviatoric thermal stresses, Sak, defined as the 
difference between the stress tensor and the nominal pressure (diagonal) tensor. In Voigt's 
notation: 



5a k 
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where H(n) is the Heaviside step function, equal to for (3 — k) strictly negative and 1 
otherwise. Deviatoric thermal stresses are caused by the vibrational (zero-point and thermal) 
energies and are shown to be important at high pressures and temperatures. The larger the 
temperature, the more visible these stresses are. 

We have previously shown that these deviatoric stresses can be relaxed to first order 
if one knows the elastic constant tensor, c^-(P, T), calculated within the (SC) QHA. 9 The 
latter are obtained from the Gibbs free energy, G, 



G(P, T) = F + PV 



(4) 



by calculating the second derivative of G with respect to the strains e. and ey^- 

,.(PT) = I^ 



(5) 



The adiabatic elastic constants, which are the relevant ones for interpretation of seismic data, 
are then computed using appropriate thermodynamics relations.- 1 ^ Below all calculated 
elastic constants, bulk and shear moduli, and velocities are adiabatic. 

Lattice parameters at high pressures and temperatures under hydrostatic conditions can 
then be corrected to first order by evaluating the strains, e&, involved in the relaxation of 
the deviatoric thermal stresses given in Eq. ([3]): 

6 

e k (P,T) = J2^(P,T)5a m . (6) 
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The Cartesian components of the relaxed lattice vectors are then: 
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are respectively the matrices of lattice vectors (a, b, c) and Cartesian strains (keeping up 
with Voigt's notation). Notice that increase in symmetry or symmetry break (phase trans- 
formations) may be induced by deviatoric thermal stresses in the presence of soft phonon, 
i.e., h and h* do not necessarily have to the same space group. 

In Ref. [9| we pointed that attainment of zero deviatoric thermal stresses within the QHA 
should involve a self- consistent cycle with simultaneous recalculation of the elastic constant 
tensor under hydrostatic condition followed by new structural relaxation, and so on. How- 
ever, such procedure is extremely computationally intensive given the need to recompute 
vibrational density of states on a PT grid every step of the cycle. We show next how to 
obtain the elastic constant tensor corrected to first order with knowledge of (EI) only. 

The components of the elastic constant tensor expanded in a Taylor series of strains (in 
Voigt's notation) defined by Eq. (jSJ) are: 
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Neglecting second and higher order terms one has: 
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In the last step above we assumed that pressure is unaffected by shear stresses, i.e., 
dP 

= for m = 4, 5, and 6, thus reducing the index summation from 6 to 3. The 

/'.:/' 

are determined using the definition of the pressure as the 
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stress derivatives of P in Eq. 



trace of the stress tensor, P = — o~ m . Taking the derivative of the pressure as function 

dP 



of each stress leads to 



da r . 
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, for m =1, 2, and 3. Therefore the first order corrected 



PT 



elastic constants at the strains given by Eq. (jo]) is reduced to 
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where 6~o~y, = ^^<5crfc. This correction requires only knowledge of the pressure derivatives of 

k=l 

Cy's which are known from the statically constrained QHA calculation, and the deviatoric 
thermal stresses given by Eq. ([3]). It gives to first order the elastic constants corrected for 
deviatoric stresses without having to explicitly calculate Gibbs free energy at the relaxed 
lattice parameters. 

As a final remark, we point that Eq. ( flOl) could also be used and tested on experimental 
data as a mean for correcting any type of deviatoric stresses, as long as the stress deviations 
remain small compared to the hydrostatic pressure (in a limit for the Taylor expansion to be 
valid). The correction only requires knowledge of (i) the three components 8<Jk, k = 1,2,3, 
and at the same time (ii) the pressure variation of the elastic constants at specified P and 
T: ^p- . Principal strain deviations, e± and e\\, are measurable quantities, for instance, 
using diffraction ring measurements 10 and their corresponding stresses are therefore also 
available from experiments. Pressure variation of the elastic constants^ are measurable 
quantities^ that require only few additional runs for estimating experimentally the pressure 
derivative of Cjj at given PT's. Eventual experimental setting that combines simultaneously 
measurements of (i) and (ii) above can be used to measure the correction to the elastic 
constants due to deviatoric stresses in DAC apparatus, after applying Eq. (fTUl) . 



III. ELASTIC CONSTANTS OF PV AND PPV 



We present in this section new results on the deviatoric thermal stresses of PPV and the 
correction to the elastic constants obtained using the (SC) QHA.- Since deviatoric thermal 
stresses of PV were recently published,- we also give here the corresponding correction to 
the elastic constants of PV. 

The PT dependent elastic constant tensors of PV and PPV determined using the (SC) 
QHA have been reported respectively in Ref. y| and Ref. [g. These are the major phases 
of the Earth's lower mantle and their elastic properties are central information for the in- 
terpretation of the seismic properties of this inaccessible region in terms of temperature, 
composition, and mineralogy. PV and PPV are both orthorhombic crystals respectively 
with symmetry Pbnm and Cmcm. This difference of symmetry group implies in partic- 
ular, as stated in Ref. [a, that "the [100]ppy, [OlOjppy, and [001]ppy directions in the 
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Cmcm structure correspond to the [110]py, [110]py, and [001]py in the Pbnm structure, 
respectively," corresponding to a rotation of 45° of the a — b reciprocal lattices. Lattice de- 
formations and deviatoric thermal stresses between PV and PPV are thus comparable only 
through this transformation. Figure (Ha) shows the deviatoric thermal stresses for PPV. 
Equivalent results for PV have recently been reported in Ref. |9| along with the analysis of 
its crystalline structure at high PT. The deviatoric stresses 5ai and 5<J2 in PPV have op- 
posite sign but similar magnitudes to that of PV (see Ref. |9|), except along the c crystalline 
axes. As stated above, deviatoric thermal stresses for PV and PPV induce distinct defor- 
mations along lattices a and b. The deviatoric thermal stresses in the z direction of PPV is 
considerably larger than the corresponding one in PV leading to larger corrections in PPV 
than in PV, as shown below. Figure QJb) shows the percentage corrections to the lattice 
parameters of PPV, based on Eq. ([6]). Interestingly, Fig. [1] shows that zero-point energy 
(the black zero Kelvin line in that figure) also produces deviatoric stresses. With increasing 
temperature, these stresses are enhanced but their origin is the anisotropic nature of the 
phonon dispersions. 

Figure [2] shows the resulting summation of the three deviatoric thermal stresses <5<j£ [of 
Fig.lJja)] for PPV (and see Ref. |9| for PV's deviatoric thermal stresses). It represents the first 
of the two ingredients necessary for the correction given by Eq. (fit)]) . Clearly, the correction 
for PPV is considerably larger than the one for PV. This is mostly due to £03 that is larger 
in PPV than in PV (see above). The correction for PPV is always negative, which has 
the effect of decreasing its elastic constants, while for PV, the correction can be negative 
(mostly at low temperature) or positive (mostly at high temperature). In principle, there 
are no reasons for having deviations of systematic nature and they should vary depending 
on the crystalline structure. One observation that remains true for all crystalline structures, 
however, is that positive deviations in one direction are to be compensated by a negative 
deviation in another direction, as observed in both PV and PPV. 

Figure [3] shows the pressure derivatives, dcij/dP, of all the elastic constants of PV and 
PPV, which is the second ingredient required for the correction according to Eq. ( FlOl) . The 
figure shows the variations of with pressure for only two temperatures, K and 3000 
K, the latter being close to the temperature of the D ' layer in the lower mantle, where the 
PPV phase is important in the geophysical models.— 

Figure H] shows the corrected bulk and shear moduli, after applying Hill's 21 (arithmetic) 
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average to the elastic constants, at several temperatures. The corrections are largest at high 
pressure and high temperature in both PV and PPV. The nature of the correction is also 
structure-dependent. Notice that the general aspect of the correction to the bulk moduli in 
Fig.H]is similar to 5a^ displayed in Fig. [21 indicating that the dominant term in the correction 
of Eq. (1101) is the deviatoric thermal stress, and to a lesser extent the pressure derivatives 
of the elastic constants. However, all corrections remain relatively small, meaning the (SC) 
QHA calculation does not suffer from significant deviatoric thermal stresses, although they 
can very well be corrected to any level of accuracy. 

Table [I] summarizes the corrections to the (SC) QHA for the elastic constants at T 
= 3000K for two pressures, P = 100 GPa and P = 120 GPa. Corrections are given in 
parenthesis. Bulk and shear moduli calculated using Voigt (uniform strain), Reuss (uniform 
stress), and Hill (arithmetic average between Voigt and Reuss) are shown. 21 The volume 
correction, abc x (1 — ei)(l — e 2 )(l — e 3 ), as shown in Fig. UJb), is reported as density, 
p(P,T). Velocities are then evaluated from Voigt-Reuss-Hill moduli, since it provides a 
realistic estimation of the true moduli. Notice that velocities are only slightly modified, 
because moduli are corrected along with the density; therefore, their ratio remains relatively 
unaltered. 

IV. CONCLUSIONS 

In summary, we have introduced a scheme to correct high PT elastic constants obtained 
using the statically constrained quasiharmonic approximation for deviatoric thermal stresses 
that develop in calculations of anisotropic structures. This self-consistent scheme was used to 
compute to first order the elastic constants of the geophysically important MgSiOs-perovskite 
and MgSi03-post-perovskite phases of the lower mantle. The corrections introduced by 
relaxation of these deviatoric stresses are quite small at relevant conditions of the lower 
mantle and previous (SC) QHA results remain essentially unchanged. However, this might 
not be the general case and the current scheme may be used to arbitrary order for computing 
high PT elastic constants to the desired level of accuracy. 
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FIG. 1: (Color online) (a) Deviatoric thermal stresses in PPV; (b) percentage lattice constant 
corrections in PPV. 5a% and 8<J2 have opposite signs and similar magnitude, similarly to the case 
of PV£ However, 5a^ in PPV is considerably larger than in PV. 
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FIG. 2: (Color online) Sum of diagonal deviatoric stresses for (a) PV and (b) PPV, as defined in 
Eq. (|10p . This sum is considerably larger in PPV because of the larger contribution from Sa^ in 
PPV. Note that the pressure ranges between PV and PPV differ, corresponding to their respective 
QHA regions of validity. 9 
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FIG. 3: (Color online) Derivatives of elastic constants as function of pressure of (a) PV and (b) 
PPV. See also note in the caption of Fig. [2j 
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TABLE I: Elastic moduli of PV— and PPV— with corrections given in parenthesis, as described 
by Eq. (jlOp . Pressure and elastic constants are in GPa, velocities in km/s, temperature in K, 
densities, /?, in g/cm 3 . The corrections are significant for bulk and shear moduli. Velocities are 
only slightly changed by the correction. Vp = \J (K H + 4/3G H )/p, Vs = \/Gh/p, V<s> = ^K H /p, 
and $ = K H / p, where upper indeces R, V, and H represent Reuss, Voigt and the Hill averages.— 
Notice that $ = V§, - 4/3 Vf. Velocit ies are calculated using the Hill averages. Decimal digits are 
presented to show the magnitude of the corrections. However, except for p, the accuracy of results 
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